AmesHousing com tidymodels

1. Introdução

Será feita a modelagem da base AmesHousing com o pacote tidymodels e alguns dos pacotes apresentados em aula ao longo do curso de Modelos Preditivos do Programa Avançado em Data Science e Decisão do Insper

Objetivo:

Modelar uma previsão do preço de venda das casas com a maior acurácia possível da base AmesHousing.

Modelos:

Linear (com ou sem seleção stepwise)
LASSO
Ridge Regression
Bagging
Floresta Aleatória

Bibliotecas

library(AmesHousing)
library(tidyverse)
library(tidymodels)
library(skimr)
library(naniar)
library(GGally)
library(vip)

Dados

Será carregada a base “make_ordinal_names”, pois que nela algumas colunas já possuem fatores ordenados, o que facilitará a modelagem. Com skim, temos uma ideia do que há na base.

dados <- make_ordinal_ames()
skim(dados)

Data summary
Name	dados
Number of rows	2930
Number of columns	81
_______________________
Column type frequency:
factor	46
numeric	35
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
MS_SubClass	0	1	FALSE	16	One: 1079, Two: 575, One: 287, One: 192
MS_Zoning	0	1	FALSE	7	Res: 2273, Res: 462, Flo: 139, Res: 27
Street	0	1	FALSE	2	Pav: 2918, Grv: 12
Alley	0	1	FALSE	3	No_: 2732, Gra: 120, Pav: 78
Lot_Shape	0	1	TRUE	4	Reg: 1859, Sli: 979, Mod: 76, Irr: 16
Land_Contour	0	1	TRUE	4	Lvl: 2633, HLS: 120, Bnk: 117, Low: 60
Utilities	0	1	TRUE	3	All: 2927, NoS: 2, NoS: 1, ELO: 0
Lot_Config	0	1	FALSE	5	Ins: 2140, Cor: 511, Cul: 180, FR2: 85
Land_Slope	0	1	TRUE	3	Gtl: 2789, Mod: 125, Sev: 16
Neighborhood	0	1	FALSE	28	Nor: 443, Col: 267, Old: 239, Edw: 194
Condition_1	0	1	FALSE	9	Nor: 2522, Fee: 164, Art: 92, RRA: 50
Condition_2	0	1	FALSE	8	Nor: 2900, Fee: 13, Art: 5, Pos: 4
Bldg_Type	0	1	FALSE	5	One: 2425, Twn: 233, Dup: 109, Twn: 101
House_Style	0	1	FALSE	8	One: 1481, Two: 873, One: 314, SLv: 128
Overall_Qual	0	1	TRUE	10	Ave: 825, Abo: 732, Goo: 602, Ver: 350
Overall_Cond	0	1	TRUE	9	Ave: 1654, Abo: 533, Goo: 390, Ver: 144
Roof_Style	0	1	FALSE	6	Gab: 2321, Hip: 551, Gam: 22, Fla: 20
Roof_Matl	0	1	FALSE	8	Com: 2887, Tar: 23, WdS: 9, WdS: 7
Exterior_1st	0	1	FALSE	16	Vin: 1026, Met: 450, HdB: 442, Wd : 420
Exterior_2nd	0	1	FALSE	17	Vin: 1015, Met: 447, HdB: 406, Wd : 397
Mas_Vnr_Type	0	1	FALSE	5	Non: 1775, Brk: 880, Sto: 249, Brk: 25
Exter_Qual	0	1	TRUE	4	Typ: 1799, Goo: 989, Exc: 107, Fai: 35
Exter_Cond	0	1	TRUE	5	Typ: 2549, Goo: 299, Fai: 67, Exc: 12
Foundation	0	1	FALSE	6	PCo: 1310, CBl: 1244, Brk: 311, Sla: 49
Bsmt_Qual	0	1	TRUE	6	Typ: 1283, Goo: 1219, Exc: 258, Fai: 88
Bsmt_Cond	0	1	TRUE	6	Typ: 2616, Goo: 122, Fai: 104, No_: 80
Bsmt_Exposure	0	1	TRUE	5	No: 1906, Av: 418, Gd: 284, Mn: 239
BsmtFin_Type_1	0	1	TRUE	7	GLQ: 859, Unf: 851, ALQ: 429, Rec: 288
BsmtFin_Type_2	0	1	TRUE	7	Unf: 2499, Rec: 106, LwQ: 89, No_: 81
Heating	0	1	FALSE	6	Gas: 2885, Gas: 27, Gra: 9, Wal: 6
Heating_QC	0	1	TRUE	5	Exc: 1495, Typ: 864, Goo: 476, Fai: 92
Central_Air	0	1	FALSE	2	Y: 2734, N: 196
Electrical	1	1	TRUE	5	SBr: 2682, Fus: 188, Fus: 50, Fus: 8
Kitchen_Qual	0	1	TRUE	5	Typ: 1494, Goo: 1160, Exc: 205, Fai: 70
Functional	0	1	TRUE	8	Typ: 2728, Min: 70, Min: 65, Mod: 35
Fireplace_Qu	0	1	TRUE	6	No_: 1422, Goo: 744, Typ: 600, Fai: 75
Garage_Type	0	1	FALSE	7	Att: 1731, Det: 782, Bui: 186, No_: 157
Garage_Finish	0	1	TRUE	4	Unf: 1231, RFn: 812, Fin: 728, No_: 159
Garage_Qual	0	1	TRUE	6	Typ: 2615, No_: 159, Fai: 124, Goo: 24
Garage_Cond	0	1	TRUE	6	Typ: 2665, No_: 159, Fai: 74, Goo: 15
Paved_Drive	0	1	TRUE	3	Pav: 2652, Dir: 216, Par: 62
Pool_QC	0	1	TRUE	5	No_: 2917, Goo: 4, Exc: 4, Typ: 3
Fence	0	1	TRUE	5	No_: 2358, Min: 330, Goo: 118, Goo: 112
Misc_Feature	0	1	FALSE	6	Non: 2824, She: 95, Gar: 5, Oth: 4
Sale_Type	0	1	FALSE	10	WD : 2536, New: 239, COD: 87, Con: 26
Sale_Condition	0	1	FALSE	6	Nor: 2413, Par: 245, Abn: 190, Fam: 46

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Lot_Frontage	1	57.65	33.50	0.00	43.00	63.00	78.00	313.00	▇▇▁▁▁
Lot_Area	1	10147.92	7880.02	1300.00	7440.25	9436.50	11555.25	215245.00	▇▁▁▁▁
Year_Built	1	1971.36	30.25	1872.00	1954.00	1973.00	2001.00	2010.00	▁▂▃▆▇
Year_Remod_Add	1	1984.27	20.86	1950.00	1965.00	1993.00	2004.00	2010.00	▅▂▂▃▇
Mas_Vnr_Area	1	101.10	178.63	0.00	0.00	0.00	162.75	1600.00	▇▁▁▁▁
BsmtFin_SF_1	1	4.18	2.23	0.00	3.00	3.00	7.00	7.00	▃▂▇▁▇
BsmtFin_SF_2	1	49.71	169.14	0.00	0.00	0.00	0.00	1526.00	▇▁▁▁▁
Bsmt_Unf_SF	1	559.07	439.54	0.00	219.00	465.50	801.75	2336.00	▇▅▂▁▁
Total_Bsmt_SF	1	1051.26	440.97	0.00	793.00	990.00	1301.50	6110.00	▇▃▁▁▁
First_Flr_SF	1	1159.56	391.89	334.00	876.25	1084.00	1384.00	5095.00	▇▃▁▁▁
Second_Flr_SF	1	335.46	428.40	0.00	0.00	0.00	703.75	2065.00	▇▃▂▁▁
Low_Qual_Fin_SF	1	4.68	46.31	0.00	0.00	0.00	0.00	1064.00	▇▁▁▁▁
Gr_Liv_Area	1	1499.69	505.51	334.00	1126.00	1442.00	1742.75	5642.00	▇▇▁▁▁
Bsmt_Full_Bath	1	0.43	0.52	0.00	0.00	0.00	1.00	3.00	▇▆▁▁▁
Bsmt_Half_Bath	1	0.06	0.25	0.00	0.00	0.00	0.00	2.00	▇▁▁▁▁
Full_Bath	1	1.57	0.55	0.00	1.00	2.00	2.00	4.00	▁▇▇▁▁
Half_Bath	1	0.38	0.50	0.00	0.00	0.00	1.00	2.00	▇▁▅▁▁
Bedroom_AbvGr	1	2.85	0.83	0.00	2.00	3.00	3.00	8.00	▁▇▂▁▁
Kitchen_AbvGr	1	1.04	0.21	0.00	1.00	1.00	1.00	3.00	▁▇▁▁▁
TotRms_AbvGrd	1	6.44	1.57	2.00	5.00	6.00	7.00	15.00	▁▇▂▁▁
Fireplaces	1	0.60	0.65	0.00	0.00	1.00	1.00	4.00	▇▇▁▁▁
Garage_Cars	1	1.77	0.76	0.00	1.00	2.00	2.00	5.00	▅▇▂▁▁
Garage_Area	1	472.66	215.19	0.00	320.00	480.00	576.00	1488.00	▃▇▃▁▁
Wood_Deck_SF	1	93.75	126.36	0.00	0.00	0.00	168.00	1424.00	▇▁▁▁▁
Open_Porch_SF	1	47.53	67.48	0.00	0.00	27.00	70.00	742.00	▇▁▁▁▁
Enclosed_Porch	1	23.01	64.14	0.00	0.00	0.00	0.00	1012.00	▇▁▁▁▁
Three_season_porch	1	2.59	25.14	0.00	0.00	0.00	0.00	508.00	▇▁▁▁▁
Screen_Porch	1	16.00	56.09	0.00	0.00	0.00	0.00	576.00	▇▁▁▁▁
Pool_Area	1	2.24	35.60	0.00	0.00	0.00	0.00	800.00	▇▁▁▁▁
Misc_Val	1	50.64	566.34	0.00	0.00	0.00	0.00	17000.00	▇▁▁▁▁
Mo_Sold	1	6.22	2.71	1.00	4.00	6.00	8.00	12.00	▅▆▇▃▃
Year_Sold	1	2007.79	1.32	2006.00	2007.00	2008.00	2009.00	2010.00	▇▇▇▇▃
Sale_Price	1	180796.06	79886.69	12789.00	129500.00	160000.00	213500.00	755000.00	▇▇▁▁▁
Longitude	1	-93.64	0.03	-93.69	-93.66	-93.64	-93.62	-93.58	▅▅▇▆▁
Latitude	1	42.03	0.02	41.99	42.02	42.03	42.05	42.06	▂▂▇▇▇

Aparentemente existe um valor faltando na coluna Electrical. Vamos carregar a base make_ames e verificar se esse valor está faltando também, novamnte com skim.

dados2 <- make_ames()
skim(dados2$Electrical)

Data summary
Name	dados2$Electrical
Number of rows	2930
Number of columns	1
_______________________
Column type frequency:
factor	1
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
data	0	1	FALSE	6	SBr: 2682, Fus: 188, Fus: 50, Fus: 8

Pelo visto não há nenhum valor faltando nessa coluna. Vamos olhar mais de perto. Primeiro encontramos a linha com valor faltando na base original, e depois procuramos o valor dessa linha na nova base.

row_na <- dados %>% 
  rowid_to_column() %>%
  filter(is.na(dados$Electrical)) %>% 
  select(rowid)

(dados$Electrical[row_na$rowid])

## [1] <NA>
## Levels: Mix < FuseP < FuseF < FuseA < SBrkr

(dados2$Electrical[row_na$rowid])

## [1] Unknown
## Levels: FuseA FuseF FuseP Mix SBrkr Unknown

O valor é “Unknown”, assim como na base ordenada. Como em ambas não há a informação de “Electrical” em apenas uma linha, vamos utilizar dropna() na base ordenada e seguir com a análise.

dados <- dados %>% 
  drop_na()

skim(dados$Electrical)

Data summary
Name	dados$Electrical
Number of rows	2929
Number of columns	1
_______________________
Column type frequency:
factor	1
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
data	0	1	TRUE	5	SBr: 2682, Fus: 188, Fus: 50, Fus: 8

2. Análise Exploratória

Para entender melhor o compartamento dos dados, será feita uma análise exploratória da base buscando relações com a variável predita “Sale_Price”.

2.1. Distribuição dos preços

Iniciamos a análise entendendo o perfil dessa variável:

Por região

dados %>%
  ggplot(aes(Longitude, Latitude, color = Sale_Price)) +
  geom_point(size = 0.5, alpha = 0.7) +
  scale_color_gradient(low = 'skyblue', high = 'darkred', labels = scales::label_number_si())

Verifica-se uma concentração de casas mais caras principalmente no norte da região analisada. Provavelmente ‘Neighborhood’ (bairro), assim como Latitude e Longitude, devem ter uma correlação alta com o preço.

Bairro

dados %>%
  ggplot(aes(Neighborhood, Sale_Price, fill = Neighborhood)) +
  geom_boxplot(show.legend = FALSE)+
  scale_y_continuous(labels = scales::label_number_si())+
  coord_flip()

Histograma

dados %>% 
    ggplot(aes(x=Sale_Price)) +
    geom_histogram(fill="skyblue", binwidth = 10000) +
    scale_x_continuous(breaks= seq(0, 800000, by=100000), labels = scales::label_number_si())

(summary(dados$Sale_Price))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   12789  129500  160000  180801  213500  755000

No histograma observa-se uma concentração mais perto do começo da amostra. Isso é esperado, pois menos pessoas conseguem comprar casas mais caras. Importante notar também que a média e mediana são próximas.

2.2. Correlação entre variáveis numéricas

Nesse item o objetivo é ter uma ideia de quais variáveis estão mais correlacionadas com “Sale_Price” para ter uma ideia do preditores e para nortear o restante da análise exploratória.

 #variáveis numéricas
num_vars <- which(sapply(dados, is.numeric))

#correlação entre todas as variáveis
corr_all <- cor(dados[,num_vars], use="pairwise.complete.obs") 

#correlação com Sale_Price em ordem decrescente
(corr_sorted <- as.matrix(sort(corr_all[,'Sale_Price'], decreasing = TRUE)))

##                           [,1]
## Sale_Price          1.00000000
## Gr_Liv_Area         0.70677666
## Garage_Cars         0.64759257
## Garage_Area         0.64013460
## Total_Bsmt_SF       0.63269326
## First_Flr_SF        0.62173389
## Year_Built          0.55861903
## Full_Bath           0.54570831
## Year_Remod_Add      0.53314636
## Mas_Vnr_Area        0.50219365
## TotRms_AbvGrd       0.49550750
## Fireplaces          0.47457710
## Wood_Deck_SF        0.32714767
## Open_Porch_SF       0.31293846
## Latitude            0.29101035
## Half_Bath           0.28520178
## Bsmt_Full_Bath      0.27580905
## Second_Flr_SF       0.26943829
## Lot_Area            0.26654763
## Lot_Frontage        0.20190876
## Bsmt_Unf_SF         0.18329078
## Bedroom_AbvGr       0.14392488
## Screen_Porch        0.11213709
## Pool_Area           0.06840003
## Mo_Sold             0.03523475
## Three_season_porch  0.03221900
## BsmtFin_SF_2        0.00600098
## Misc_Val           -0.01569664
## Year_Sold          -0.03056032
## Bsmt_Half_Bath     -0.03583132
## Low_Qual_Fin_SF    -0.03766575
## Kitchen_AbvGr      -0.11982695
## Enclosed_Porch     -0.12881128
## BsmtFin_SF_1       -0.13487107
## Longitude          -0.25141580

#correlação alta com Sale_Price
corr_high <- names(which(apply(corr_sorted, 1, function(x) abs(x)>0.5))) 

dados[,corr_high] %>% 
  ggpairs()

2.3. Ano de construção, Área total e Material de construção

dados %>%
  ggplot(aes(Year_Built, Sale_Price, size = Lot_Area/1000,  color = Foundation)) +
  geom_point(alpha = 0.5)+
  scale_y_continuous(labels = scales::label_number_si())

O mais marcante nesse gráfico é a mudança do material utilizado na estrutura das casas ao longo dos anos. A relação entre ano de construção e preço é também confirmada, visto que ‘Year_Built’ é umas das variáveis numéricas com correlação maior do que 0.5.

2.4. Garagem

No topo da lista de variáveis correlacionadas está a a área da garagem e quantidade de carros. Vamos observar o comportamento destas e outras variáveis relacionadas.

dados %>%
  ggplot(aes(Garage_Area, Sale_Price, size = Garage_Cars, color = Garage_Type)) +
  geom_point(alpha = 0.5)+
  scale_y_continuous(labels = scales::label_number_si())

dados %>%
  ggplot(aes(Garage_Qual, Sale_Price, color = Garage_Finish)) +
  geom_boxplot()+
  scale_y_continuous(labels = scales::label_number_si())

Além das variáveis altamente correlacioanadas, há uma concentração do tipo de garagem em difenrentes níveis de preço. O tipo de acabamento, no entanto, parece não ter nenhuma relação significativa.

2.5. Área total, Lareiras e Condição de Venda

dados %>%
  ggplot(aes(Gr_Liv_Area, Sale_Price, size = Fireplaces, color = Sale_Condition)) +
  geom_point(alpha = 0.5)+
  scale_y_continuous(labels = scales::label_number_si())

Aqui apenas é confirmada a correlação da quantidade de lareiras e área comum. A condição de venda “Normal” predomina na amostra.

2.5 Qualidade de venda

dados %>%
  ggplot(aes(Overall_Qual, Sale_Price, fill = Street)) +
  geom_boxplot()+
  scale_y_continuous(labels = scales::label_number_si())

Apesar de estar na base como uma variável categórica, “Overall Quality” parace ter uma correlação alta com Sale_Price.

3. Pré processamento

Antes de alimentar os dados aos modelos, é necessário separar a base em treino e teste, e garantir que está em um formato interpretável para os algoritmos utilizados. O parâmetro “strata” divide a base proporcionalmente conforme a variável indicada.

3.1 Separação da base em treino e teste (treino/teste/total):

set.seed(123)
(ames_split <- initial_split(dados, prop = 0.8, strata = 'Sale_Price'))

## <2345/584/2929>

ames_train <- training(ames_split)
ames_test <- testing(ames_split)

3.2. Receita

A interface do tidymodels permite a criação de uma receita com o pacote “recipes” para processar os dados antes do modelo. Isso facilita também o pré-processamento da base de teste ao final do relatório, pois aplica exatamente as mesmas modificações. Algumas variáveis já estão na base como fator, portanto não é necessário ordená-las. No entanto, é necessário incluir na receita um passo para converter a ordem em números: step_ordinalscore

A seguir as variáveis ordenadas:

ord_vars <- vapply(dados, is.ordered, logical(1))

(ordered <-names(ord_vars)[ord_vars])

##  [1] "Lot_Shape"      "Land_Contour"   "Utilities"      "Land_Slope"    
##  [5] "Overall_Qual"   "Overall_Cond"   "Exter_Qual"     "Exter_Cond"    
##  [9] "Bsmt_Qual"      "Bsmt_Cond"      "Bsmt_Exposure"  "BsmtFin_Type_1"
## [13] "BsmtFin_Type_2" "Heating_QC"     "Electrical"     "Kitchen_Qual"  
## [17] "Functional"     "Fireplace_Qu"   "Garage_Finish"  "Garage_Qual"   
## [21] "Garage_Cond"    "Paved_Drive"    "Pool_QC"        "Fence"

Como algumas colunas possuem muitas categorias, as que possuem mais de 10 serão reduzidas com “step_other” (10 é o número de classificações atribuídas nas colunas que avaliam a qualidade). Após isso, cria-se variáveis dummy para as demais categorias (step_dummy), remove-se as variáveis com apenas um valor (step_vz) e normaliza-se a base (step_normalize).

(ames_rec <- recipe(Sale_Price ~ ., data = ames_train) %>%
  #step_log(Sale_Price)) %>% 
  step_other(MS_SubClass, Neighborhood, Exterior_1st, Exterior_2nd, threshold = 0.02) %>% 
  step_ordinalscore(ordered) %>%
  step_dummy(all_nominal(), -all_outcomes()) %>%
  step_zv(all_predictors()) %>% 
  step_normalize(all_numeric()) %>% 
  prep())

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         80
## 
## Training data contained 2345 data points and no missing data.
## 
## Operations:
## 
## Collapsing factor levels for MS_SubClass, Neighborhood, ... [trained]
## Scoring for Lot_Shape, Land_Contour, Utilities, ... [trained]
## Dummy variables from MS_SubClass, MS_Zoning, Street, Alley, ... [trained]
## Zero variance filter removed Roof_Matl_Membran [trained]
## Centering and scaling for Lot_Frontage, Lot_Area, Lot_Shape, ... [trained]

A receita pronta é então aplicada nas bases de treino e teste com a função “bake”. Para a base teste, pode-se utilizar juice, que é um caso especial de bake.

train_baked  <- juice(ames_rec) 

test_baked <- bake(ames_rec, new_data = ames_test)

A seguir uma ideia de como ficou a base processada:

skim(train_baked)

Data summary
Name	train_baked
Number of rows	2345
Number of columns	190
_______________________
Column type frequency:
numeric	190
________________________
Group variables	None

Variable type: numeric

skim_variable	complete_rate	sd	p0	p25	p50	p75	p100	hist
Lot_Frontage	1	1	-1.75	-0.42	0.15	0.64	7.73	▇▇▁▁▁
Lot_Area	1	1	-1.10	-0.33	-0.08	0.19	25.71	▇▁▁▁▁
Lot_Shape	1	1	-4.61	-1.07	0.70	0.70	0.70	▁▁▁▅▇
Land_Contour	1	1	-4.67	0.31	0.31	0.31	0.31	▁▁▁▁▇
Utilities	1	1	-34.22	0.03	0.03	0.03	0.03	▁▁▁▁▇
Land_Slope	1	1	-7.87	0.22	0.22	0.22	0.22	▁▁▁▁▇
Overall_Qual	1	1	-3.62	-0.77	-0.06	0.65	2.79	▁▂▇▅▁
Overall_Cond	1	1	-4.18	-0.52	-0.52	0.40	3.15	▁▁▇▅▁
Year_Built	1	1	-3.25	-0.59	0.06	0.95	1.28	▁▂▃▆▇
Year_Remod_Add	1	1	-1.64	-0.92	0.43	0.91	1.24	▅▂▂▃▇
Mas_Vnr_Area	1	1	-0.57	-0.57	-0.57	0.36	6.84	▇▁▁▁▁
Exter_Qual	1	1	-2.40	-0.68	-0.68	1.05	2.77	▁▇▁▅▁
Exter_Cond	1	1	-5.72	-0.23	-0.23	-0.23	5.26	▁▁▇▁▁
Bsmt_Qual	1	1	-3.87	-0.53	0.58	0.58	1.69	▁▁▇▇▂
Bsmt_Cond	1	1	-5.14	0.12	0.12	0.12	3.63	▁▁▇▁▁
Bsmt_Exposure	1	1	-1.52	-0.58	-0.58	0.35	2.23	▁▇▁▂▁
BsmtFin_Type_1	1	1	-1.69	-1.22	0.20	1.15	1.15	▆▁▂▂▇
BsmtFin_SF_1	1	1	-1.41	-0.52	-0.52	1.27	1.27	▅▆▁▁▇
BsmtFin_Type_2	1	1	-1.33	-0.30	-0.30	-0.30	4.86	▇▁▁▁▁
BsmtFin_SF_2	1	1	-0.30	-0.30	-0.30	-0.30	8.58	▇▁▁▁▁
Bsmt_Unf_SF	1	1	-1.27	-0.77	-0.21	0.55	4.07	▇▅▂▁▁
Total_Bsmt_SF	1	1	-2.46	-0.61	-0.14	0.58	9.47	▇▇▁▁▁
Heating_QC	1	1	-3.29	-1.20	0.89	0.89	0.89	▁▁▅▂▇
Electrical	1	1	-9.89	0.28	0.28	0.28	0.28	▁▁▁▁▇
First_Flr_SF	1	1	-2.13	-0.73	-0.19	0.59	10.22	▇▃▁▁▁
Second_Flr_SF	1	1	-0.79	-0.79	-0.79	0.85	4.04	▇▃▂▁▁
Low_Qual_Fin_SF	1	1	-0.10	-0.10	-0.10	-0.10	22.30	▇▁▁▁▁
Gr_Liv_Area	1	1	-2.35	-0.75	-0.12	0.49	7.26	▅▇▁▁▁
Bsmt_Full_Bath	1	1	-0.82	-0.82	-0.82	1.08	4.88	▇▆▁▁▁
Bsmt_Half_Bath	1	1	-0.25	-0.25	-0.25	-0.25	7.84	▇▁▁▁▁
Full_Bath	1	1	-2.81	-1.01	0.78	0.78	4.38	▁▇▇▁▁
Half_Bath	1	1	-0.76	-0.76	-0.76	1.22	3.19	▇▁▅▁▁
Bedroom_AbvGr	1	1	-3.41	-1.02	0.18	0.18	6.16	▁▇▂▁▁
Kitchen_AbvGr	1	1	-4.92	-0.20	-0.20	-0.20	9.22	▁▇▁▁▁
Kitchen_Qual	1	1	-2.28	-0.76	-0.76	0.76	2.28	▁▇▁▆▁
TotRms_AbvGrd	1	1	-2.80	-0.91	-0.28	0.35	5.39	▁▇▂▁▁
Functional	1	1	-10.53	0.24	0.24	0.24	0.24	▁▁▁▁▇
Fireplaces	1	1	-0.93	-0.93	0.61	0.61	5.23	▇▇▁▁▁
Fireplace_Qu	1	1	-0.99	-0.99	0.11	1.22	1.77	▇▁▃▅▁
Garage_Finish	1	1	-1.91	-0.79	0.32	0.32	1.44	▁▇▁▆▅
Garage_Cars	1	1	-2.31	-1.00	0.31	0.31	4.23	▁▃▇▂▁
Garage_Area	1	1	-2.20	-0.70	0.04	0.49	4.75	▃▇▂▁▁
Garage_Qual	1	1	-3.87	0.28	0.28	0.28	3.05	▁▁▁▇▁
Garage_Cond	1	1	-3.90	0.27	0.27	0.27	3.05	▁▁▁▇▁
Paved_Drive	1	1	-3.50	0.31	0.31	0.31	0.31	▁▁▁▁▇
Wood_Deck_SF	1	1	-0.75	-0.75	-0.75	0.61	6.26	▇▂▁▁▁
Open_Porch_SF	1	1	-0.70	-0.70	-0.32	0.33	10.20	▇▁▁▁▁
Enclosed_Porch	1	1	-0.38	-0.38	-0.38	-0.38	8.72	▇▁▁▁▁
Three_season_porch	1	1	-0.10	-0.10	-0.10	-0.10	19.11	▇▁▁▁▁
Screen_Porch	1	1	-0.29	-0.29	-0.29	-0.29	9.92	▇▁▁▁▁
Pool_Area	1	1	-0.06	-0.06	-0.06	-0.06	21.31	▇▁▁▁▁
Pool_QC	1	1	-0.07	-0.07	-0.07	-0.07	18.33	▇▁▁▁▁
Fence	1	1	-0.48	-0.48	-0.48	-0.48	2.76	▇▁▁▁▁
Misc_Val	1	1	-0.09	-0.09	-0.09	-0.09	27.08	▇▁▁▁▁
Mo_Sold	1	1	-1.92	-0.82	-0.08	0.65	2.13	▃▆▇▃▃
Year_Sold	1	1	-1.36	-0.59	0.17	0.93	1.69	▇▇▇▇▃
Longitude	1	1	-1.94	-0.67	0.04	0.80	2.56	▅▅▇▆▁
Latitude	1	1	-2.63	-0.68	0.01	0.83	1.58	▂▂▇▇▇
Sale_Price	1	1	-2.09	-0.64	-0.26	0.40	7.14	▇▇▁▁▁
MS_SubClass_One_Story_1945_and_Older	1	1	-0.22	-0.22	-0.22	-0.22	4.55	▇▁▁▁▁
MS_SubClass_One_and_Half_Story_Finished_All_Ages	1	1	-0.34	-0.34	-0.34	-0.34	2.97	▇▁▁▁▁
MS_SubClass_Two_Story_1946_and_Newer	1	1	-0.50	-0.50	-0.50	-0.50	2.01	▇▁▁▁▂
MS_SubClass_Two_Story_1945_and_Older	1	1	-0.22	-0.22	-0.22	-0.22	4.57	▇▁▁▁▁
MS_SubClass_Split_or_Multilevel	1	1	-0.20	-0.20	-0.20	-0.20	5.09	▇▁▁▁▁
MS_SubClass_Duplex_All_Styles_and_Ages	1	1	-0.19	-0.19	-0.19	-0.19	5.19	▇▁▁▁▁
MS_SubClass_One_Story_PUD_1946_and_Newer	1	1	-0.27	-0.27	-0.27	-0.27	3.75	▇▁▁▁▁
MS_SubClass_Two_Story_PUD_1946_and_Newer	1	1	-0.21	-0.21	-0.21	-0.21	4.66	▇▁▁▁▁
MS_SubClass_Two_Family_conversion_All_Styles_and_Ages	1	1	-0.15	-0.15	-0.15	-0.15	6.84	▇▁▁▁▁
MS_SubClass_other	1	1	-0.20	-0.20	-0.20	-0.20	4.95	▇▁▁▁▁
MS_Zoning_Residential_High_Density	1	1	-0.10	-0.10	-0.10	-0.10	10.27	▇▁▁▁▁
MS_Zoning_Residential_Low_Density	1	1	-1.85	0.54	0.54	0.54	0.54	▂▁▁▁▇
MS_Zoning_Residential_Medium_Density	1	1	-0.43	-0.43	-0.43	-0.43	2.31	▇▁▁▁▂
MS_Zoning_A_agr	1	1	-0.03	-0.03	-0.03	-0.03	34.22	▇▁▁▁▁
MS_Zoning_C_all	1	1	-0.10	-0.10	-0.10	-0.10	10.27	▇▁▁▁▁
MS_Zoning_I_all	1	1	-0.03	-0.03	-0.03	-0.03	34.22	▇▁▁▁▁
Street_Pave	1	1	-15.28	0.07	0.07	0.07	0.07	▁▁▁▁▇
Alley_No_Alley_Access	1	1	-3.65	0.27	0.27	0.27	0.27	▁▁▁▁▇
Alley_Paved	1	1	-0.16	-0.16	-0.16	-0.16	6.12	▇▁▁▁▁
Lot_Config_CulDSac	1	1	-0.26	-0.26	-0.26	-0.26	3.89	▇▁▁▁▁
Lot_Config_FR2	1	1	-0.17	-0.17	-0.17	-0.17	5.83	▇▁▁▁▁
Lot_Config_FR3	1	1	-0.07	-0.07	-0.07	-0.07	15.28	▇▁▁▁▁
Lot_Config_Inside	1	1	-1.65	-1.65	0.61	0.61	0.61	▃▁▁▁▇
Neighborhood_College_Creek	1	1	-0.32	-0.32	-0.32	-0.32	3.15	▇▁▁▁▁
Neighborhood_Old_Town	1	1	-0.30	-0.30	-0.30	-0.30	3.37	▇▁▁▁▁
Neighborhood_Edwards	1	1	-0.27	-0.27	-0.27	-0.27	3.71	▇▁▁▁▁
Neighborhood_Somerset	1	1	-0.26	-0.26	-0.26	-0.26	3.82	▇▁▁▁▁
Neighborhood_Northridge_Heights	1	1	-0.24	-0.24	-0.24	-0.24	4.11	▇▁▁▁▁
Neighborhood_Gilbert	1	1	-0.24	-0.24	-0.24	-0.24	4.18	▇▁▁▁▁
Neighborhood_Sawyer	1	1	-0.23	-0.23	-0.23	-0.23	4.38	▇▁▁▁▁
Neighborhood_Northwest_Ames	1	1	-0.22	-0.22	-0.22	-0.22	4.57	▇▁▁▁▁
Neighborhood_Sawyer_West	1	1	-0.22	-0.22	-0.22	-0.22	4.59	▇▁▁▁▁
Neighborhood_Mitchell	1	1	-0.19	-0.19	-0.19	-0.19	5.19	▇▁▁▁▁
Neighborhood_Brookside	1	1	-0.20	-0.20	-0.20	-0.20	4.89	▇▁▁▁▁
Neighborhood_Crawford	1	1	-0.19	-0.19	-0.19	-0.19	5.22	▇▁▁▁▁
Neighborhood_Iowa_DOT_and_Rail_Road	1	1	-0.18	-0.18	-0.18	-0.18	5.54	▇▁▁▁▁
Neighborhood_Timberland	1	1	-0.16	-0.16	-0.16	-0.16	6.39	▇▁▁▁▁
Neighborhood_Northridge	1	1	-0.15	-0.15	-0.15	-0.15	6.45	▇▁▁▁▁
Neighborhood_other	1	1	-0.34	-0.34	-0.34	-0.34	2.90	▇▁▁▁▁
Condition_1_Feedr	1	1	-0.25	-0.25	-0.25	-0.25	4.03	▇▁▁▁▁
Condition_1_Norm	1	1	-2.49	0.40	0.40	0.40	0.40	▁▁▁▁▇
Condition_1_PosA	1	1	-0.08	-0.08	-0.08	-0.08	12.46	▇▁▁▁▁
Condition_1_PosN	1	1	-0.11	-0.11	-0.11	-0.11	8.78	▇▁▁▁▁
Condition_1_RRAe	1	1	-0.09	-0.09	-0.09	-0.09	11.06	▇▁▁▁▁
Condition_1_RRAn	1	1	-0.14	-0.14	-0.14	-0.14	7.40	▇▁▁▁▁
Condition_1_RRNe	1	1	-0.05	-0.05	-0.05	-0.05	21.63	▇▁▁▁▁
Condition_1_RRNn	1	1	-0.05	-0.05	-0.05	-0.05	19.74	▇▁▁▁▁
Condition_2_Feedr	1	1	-0.07	-0.07	-0.07	-0.07	15.28	▇▁▁▁▁
Condition_2_Norm	1	1	-10.78	0.09	0.09	0.09	0.09	▁▁▁▁▇
Condition_2_PosA	1	1	-0.04	-0.04	-0.04	-0.04	27.93	▇▁▁▁▁
Condition_2_PosN	1	1	-0.03	-0.03	-0.03	-0.03	34.22	▇▁▁▁▁
Condition_2_RRAe	1	1	-0.02	-0.02	-0.02	-0.02	48.40	▇▁▁▁▁
Condition_2_RRAn	1	1	-0.02	-0.02	-0.02	-0.02	48.40	▇▁▁▁▁
Condition_2_RRNn	1	1	-0.02	-0.02	-0.02	-0.02	48.40	▇▁▁▁▁
Bldg_Type_TwoFmCon	1	1	-0.15	-0.15	-0.15	-0.15	6.71	▇▁▁▁▁
Bldg_Type_Duplex	1	1	-0.19	-0.19	-0.19	-0.19	5.19	▇▁▁▁▁
Bldg_Type_Twnhs	1	1	-0.19	-0.19	-0.19	-0.19	5.29	▇▁▁▁▁
Bldg_Type_TwnhsE	1	1	-0.30	-0.30	-0.30	-0.30	3.37	▇▁▁▁▁
House_Style_One_and_Half_Unf	1	1	-0.09	-0.09	-0.09	-0.09	11.70	▇▁▁▁▁
House_Style_One_Story	1	1	-1.00	-1.00	1.00	1.00	1.00	▇▁▁▁▇
House_Style_SFoyer	1	1	-0.17	-0.17	-0.17	-0.17	5.83	▇▁▁▁▁
House_Style_SLvl	1	1	-0.20	-0.20	-0.20	-0.20	4.89	▇▁▁▁▁
House_Style_Two_and_Half_Fin	1	1	-0.05	-0.05	-0.05	-0.05	19.74	▇▁▁▁▁
House_Style_Two_and_Half_Unf	1	1	-0.09	-0.09	-0.09	-0.09	10.78	▇▁▁▁▁
House_Style_Two_Story	1	1	-0.65	-0.65	-0.65	1.53	1.53	▇▁▁▁▃
Roof_Style_Gable	1	1	-1.95	0.51	0.51	0.51	0.51	▂▁▁▁▇
Roof_Style_Gambrel	1	1	-0.09	-0.09	-0.09	-0.09	11.70	▇▁▁▁▁
Roof_Style_Hip	1	1	-0.48	-0.48	-0.48	-0.48	2.08	▇▁▁▁▂
Roof_Style_Mansard	1	1	-0.06	-0.06	-0.06	-0.06	16.11	▇▁▁▁▁
Roof_Style_Shed	1	1	-0.05	-0.05	-0.05	-0.05	21.63	▇▁▁▁▁
Roof_Matl_CompShg	1	1	-8.78	0.11	0.11	0.11	0.11	▁▁▁▁▇
Roof_Matl_Metal	1	1	-0.02	-0.02	-0.02	-0.02	48.40	▇▁▁▁▁
Roof_Matl_Roll	1	1	-0.02	-0.02	-0.02	-0.02	48.40	▇▁▁▁▁
Roof_Matl_Tar.Grv	1	1	-0.09	-0.09	-0.09	-0.09	10.78	▇▁▁▁▁
Roof_Matl_WdShake	1	1	-0.04	-0.04	-0.04	-0.04	24.19	▇▁▁▁▁
Roof_Matl_WdShngl	1	1	-0.04	-0.04	-0.04	-0.04	24.19	▇▁▁▁▁
Exterior_1st_CemntBd	1	1	-0.22	-0.22	-0.22	-0.22	4.62	▇▁▁▁▁
Exterior_1st_HdBoard	1	1	-0.41	-0.41	-0.41	-0.41	2.44	▇▁▁▁▂
Exterior_1st_MetalSd	1	1	-0.43	-0.43	-0.43	-0.43	2.34	▇▁▁▁▂
Exterior_1st_Plywood	1	1	-0.29	-0.29	-0.29	-0.29	3.48	▇▁▁▁▁
Exterior_1st_VinylSd	1	1	-0.74	-0.74	-0.74	1.36	1.36	▇▁▁▁▅
Exterior_1st_Wd.Sdng	1	1	-0.41	-0.41	-0.41	-0.41	2.42	▇▁▁▁▂
Exterior_1st_WdShing	1	1	-0.14	-0.14	-0.14	-0.14	6.92	▇▁▁▁▁
Exterior_1st_other	1	1	-0.19	-0.19	-0.19	-0.19	5.25	▇▁▁▁▁
Exterior_2nd_HdBoard	1	1	-0.39	-0.39	-0.39	-0.39	2.59	▇▁▁▁▁
Exterior_2nd_MetalSd	1	1	-0.43	-0.43	-0.43	-0.43	2.34	▇▁▁▁▂
Exterior_2nd_Plywood	1	1	-0.32	-0.32	-0.32	-0.32	3.10	▇▁▁▁▁
Exterior_2nd_VinylSd	1	1	-0.73	-0.73	-0.73	1.37	1.37	▇▁▁▁▅
Exterior_2nd_Wd.Sdng	1	1	-0.40	-0.40	-0.40	-0.40	2.49	▇▁▁▁▁
Exterior_2nd_Wd.Shng	1	1	-0.17	-0.17	-0.17	-0.17	5.87	▇▁▁▁▁
Exterior_2nd_other	1	1	-0.26	-0.26	-0.26	-0.26	3.89	▇▁▁▁▁
Mas_Vnr_Type_BrkFace	1	1	-0.64	-0.64	-0.64	1.55	1.55	▇▁▁▁▃
Mas_Vnr_Type_CBlock	1	1	-0.02	-0.02	-0.02	-0.02	48.40	▇▁▁▁▁
Mas_Vnr_Type_None	1	1	-1.25	-1.25	0.80	0.80	0.80	▅▁▁▁▇
Mas_Vnr_Type_Stone	1	1	-0.31	-0.31	-0.31	-0.31	3.25	▇▁▁▁▁
Foundation_CBlock	1	1	-0.86	-0.86	-0.86	1.17	1.17	▇▁▁▁▆
Foundation_PConc	1	1	-0.90	-0.90	-0.90	1.11	1.11	▇▁▁▁▆
Foundation_Slab	1	1	-0.12	-0.12	-0.12	-0.12	8.12	▇▁▁▁▁
Foundation_Stone	1	1	-0.05	-0.05	-0.05	-0.05	18.27	▇▁▁▁▁
Foundation_Wood	1	1	-0.04	-0.04	-0.04	-0.04	24.19	▇▁▁▁▁
Heating_GasA	1	1	-7.69	0.13	0.13	0.13	0.13	▁▁▁▁▇
Heating_GasW	1	1	-0.10	-0.10	-0.10	-0.10	9.63	▇▁▁▁▁
Heating_Grav	1	1	-0.05	-0.05	-0.05	-0.05	18.27	▇▁▁▁▁
Heating_OthW	1	1	-0.03	-0.03	-0.03	-0.03	34.22	▇▁▁▁▁
Heating_Wall	1	1	-0.04	-0.04	-0.04	-0.04	24.19	▇▁▁▁▁
Central_Air_Y	1	1	-3.62	0.28	0.28	0.28	0.28	▁▁▁▁▇
Garage_Type_Basment	1	1	-0.11	-0.11	-0.11	-0.11	8.93	▇▁▁▁▁
Garage_Type_BuiltIn	1	1	-0.26	-0.26	-0.26	-0.26	3.78	▇▁▁▁▁
Garage_Type_CarPort	1	1	-0.07	-0.07	-0.07	-0.07	14.56	▇▁▁▁▁
Garage_Type_Detchd	1	1	-0.60	-0.60	-0.60	1.67	1.67	▇▁▁▁▃
Garage_Type_More_Than_Two_Types	1	1	-0.09	-0.09	-0.09	-0.09	10.78	▇▁▁▁▁
Garage_Type_No_Garage	1	1	-0.24	-0.24	-0.24	-0.24	4.14	▇▁▁▁▁
Misc_Feature_Gar2	1	1	-0.05	-0.05	-0.05	-0.05	21.63	▇▁▁▁▁
Misc_Feature_None	1	1	-5.09	0.20	0.20	0.20	0.20	▁▁▁▁▇
Misc_Feature_Othr	1	1	-0.04	-0.04	-0.04	-0.04	27.93	▇▁▁▁▁
Misc_Feature_Shed	1	1	-0.18	-0.18	-0.18	-0.18	5.43	▇▁▁▁▁
Misc_Feature_TenC	1	1	-0.02	-0.02	-0.02	-0.02	48.40	▇▁▁▁▁
Sale_Type_Con	1	1	-0.04	-0.04	-0.04	-0.04	24.19	▇▁▁▁▁
Sale_Type_ConLD	1	1	-0.10	-0.10	-0.10	-0.10	10.27	▇▁▁▁▁
Sale_Type_ConLI	1	1	-0.05	-0.05	-0.05	-0.05	18.27	▇▁▁▁▁
Sale_Type_ConLw	1	1	-0.05	-0.05	-0.05	-0.05	18.27	▇▁▁▁▁
Sale_Type_CWD	1	1	-0.07	-0.07	-0.07	-0.07	13.94	▇▁▁▁▁
Sale_Type_New	1	1	-0.30	-0.30	-0.30	-0.30	3.35	▇▁▁▁▁
Sale_Type_Oth	1	1	-0.05	-0.05	-0.05	-0.05	19.74	▇▁▁▁▁
Sale_Type_VWD	1	1	-0.02	-0.02	-0.02	-0.02	48.40	▇▁▁▁▁
Sale_Type_WD.	1	1	-2.51	0.40	0.40	0.40	0.40	▁▁▁▁▇
Sale_Condition_AdjLand	1	1	-0.07	-0.07	-0.07	-0.07	13.94	▇▁▁▁▁
Sale_Condition_Alloca	1	1	-0.09	-0.09	-0.09	-0.09	11.70	▇▁▁▁▁
Sale_Condition_Family	1	1	-0.13	-0.13	-0.13	-0.13	7.69	▇▁▁▁▁
Sale_Condition_Normal	1	1	-2.13	0.47	0.47	0.47	0.47	▂▁▁▁▇
Sale_Condition_Partial	1	1	-0.30	-0.30	-0.30	-0.30	3.30	▇▁▁▁▁

Aplica-se então cross validation na base de treino. Além disso, serão criadas amostras bootstrap para ajuste dos hiperparâmetros de alguns dos modelos utilizados.

set.seed(123)
(cv_splits <- vfold_cv(train_baked, v = 5, strata = Sale_Price))

## #  5-fold cross-validation using stratification 
## # A tibble: 5 x 2
##   splits             id   
##   <named list>       <chr>
## 1 <split [1.9K/470]> Fold1
## 2 <split [1.9K/470]> Fold2
## 3 <split [1.9K/469]> Fold3
## 4 <split [1.9K/468]> Fold4
## 5 <split [1.9K/468]> Fold5

set.seed(123)
(ames_boot <- bootstraps(train_baked, times = 10, strata = Sale_Price))

## # Bootstrap sampling using stratification 
## # A tibble: 10 x 2
##    splits             id         
##    <named list>       <chr>      
##  1 <split [2.3K/865]> Bootstrap01
##  2 <split [2.3K/857]> Bootstrap02
##  3 <split [2.3K/890]> Bootstrap03
##  4 <split [2.3K/858]> Bootstrap04
##  5 <split [2.3K/872]> Bootstrap05
##  6 <split [2.3K/864]> Bootstrap06
##  7 <split [2.3K/855]> Bootstrap07
##  8 <split [2.3K/856]> Bootstrap08
##  9 <split [2.3K/844]> Bootstrap09
## 10 <split [2.3K/876]> Bootstrap10

4. Modelos

Antes de ajustar os modelos, é necessário especificar o pacote usado, modo e parâmetros. Depois disso, é feita a modelagem nos diferentes folds criados por cross validation para avaliar seu desempenho.

4.1. Regressão Linear

4.1.1. Sem seleção stepwise

lm_spec <- linear_reg() %>% 
  set_engine('lm')

#fit nos cv folds
lm_res <- fit_resamples(Sale_Price ~ .,
               lm_spec,
               cv_splits,
               control = control_resamples(save_pred = TRUE))

Resumo do modelo ajustado:

(lm_res %>% 
  collect_metrics())

## # A tibble: 2 x 5
##   .metric .estimator  mean     n std_err
##   <chr>   <chr>      <dbl> <int>   <dbl>
## 1 rmse    standard   0.335     5  0.0336
## 2 rsq     standard   0.887     5  0.0197

Modelo ajustado na base completa, resumo e importância das variáveis:

lm_fit <- lm_spec %>% 
      fit(Sale_Price ~.,
      data = train_baked)

lm_fit %>% 
  summary()

##         Length Class      Mode   
## lvl      0     -none-     NULL   
## spec     5     linear_reg list   
## fit     12     lm         list   
## preproc  1     -none-     list   
## elapsed  5     proc_time  numeric

lm_fit %>% 
  tidy()

## # A tibble: 187 x 5
##    term          estimate std.error statistic  p.value
##    <chr>            <dbl>     <dbl>     <dbl>    <dbl>
##  1 (Intercept)  -1.28e-14   0.00580 -2.21e-12 1.00e+ 0
##  2 Lot_Frontage  1.95e- 2   0.00720  2.71e+ 0 6.74e- 3
##  3 Lot_Area      4.44e- 2   0.00776  5.72e+ 0 1.20e- 8
##  4 Lot_Shape    -5.80e- 3   0.00738 -7.86e- 1 4.32e- 1
##  5 Land_Contour -8.18e- 3   0.00803 -1.02e+ 0 3.08e- 1
##  6 Utilities     1.75e- 3   0.00712  2.47e- 1 8.05e- 1
##  7 Land_Slope    2.35e- 3   0.00829  2.84e- 1 7.77e- 1
##  8 Overall_Qual  1.48e- 1   0.0128   1.15e+ 1 7.32e-30
##  9 Overall_Cond  7.26e- 2   0.00852  8.52e+ 0 3.07e-17
## 10 Year_Built    1.54e- 1   0.0221   6.98e+ 0 4.05e-12
## # … with 177 more rows

4.1.2 Com seleção stepwise backward

O pacote tidymodels não possui uma integração pronta com o pacote leaps para realizar stepwise selection. Devido a isso, o procedimento será feito, mas o modelo não será considerado na avaliação final. É possível criar essa integração pelo pacote parsnip.

library(leaps)

stepb <- regsubsets(Sale_Price ~ ., data = train_baked, nvmax = 10, 
                          method = "backward")

## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax, force.in =
## force.in, : 3 linear dependencies found

## Reordering variables and trying again:

  resumo_stepb <- summary(stepb)
  #resumo_stepb$outmat
  
  coef(stepb, id = which.min(resumo_stepb$bic))

##       (Intercept)      Overall_Qual        Exter_Qual       Bsmt_Unf_SF 
##     -2.632578e-17      3.166017e-01      2.124920e-01     -1.651529e-01 
##     Total_Bsmt_SF      First_Flr_SF     Second_Flr_SF Misc_Feature_Shed 
##      2.491390e-01      3.179770e-01      3.014907e-01     -3.175356e-03 
## Misc_Feature_TenC     Sale_Type_Con   Sale_Type_ConLD     Sale_Type_WD. 
##     -9.898401e-03      3.379658e-03     -1.065269e-02     -5.593790e-02

  plot(stepb, scale = "adjr2")

4.2. Lasso

Nesse modelo, é utilizada a função tune() para encontrar o melhor lambda (penalty) nas amostras bootstrap.

lasso_spec <- linear_reg(penalty = tune(), mixture = 1) %>%
  set_engine("glmnet", standardize = FALSE) #false, pois a base já foi pré-processada anteriormente 

lambda_grid <- grid_regular(penalty(), levels = 500)

doParallel::registerDoParallel() #processamento em paralelo para otimizar

set.seed(123)
lasso_grid <- tune_grid(Sale_Price ~.,
  model = lasso_spec,
  resamples = ames_boot,
  grid = lambda_grid)

Avaliação dos parâmetros - lambda para o menor erro quadrado médio (rmse):

lasso_grid %>%
  collect_metrics()

## # A tibble: 1,000 x 6
##     penalty .metric .estimator  mean     n std_err
##       <dbl> <chr>   <chr>      <dbl> <int>   <dbl>
##  1 1.00e-10 rmse    standard   0.370    10 0.0128 
##  2 1.00e-10 rsq     standard   0.870    10 0.00796
##  3 1.05e-10 rmse    standard   0.370    10 0.0128 
##  4 1.05e-10 rsq     standard   0.870    10 0.00796
##  5 1.10e-10 rmse    standard   0.370    10 0.0128 
##  6 1.10e-10 rsq     standard   0.870    10 0.00796
##  7 1.15e-10 rmse    standard   0.370    10 0.0128 
##  8 1.15e-10 rsq     standard   0.870    10 0.00796
##  9 1.20e-10 rmse    standard   0.370    10 0.0128 
## 10 1.20e-10 rsq     standard   0.870    10 0.00796
## # … with 990 more rows

lasso_grid %>%
  collect_metrics() %>%
  ggplot(aes(penalty, mean, color = .metric)) +
  geom_errorbar(aes(
    ymin = mean - std_err,
    ymax = mean + std_err
    ),
  alpha = 0.5
  ) +
  geom_line(size = 1.5) +
  facet_wrap(~.metric, scales = "free", nrow = 2) +
  scale_x_log10() +
  labs(x = 'Log(lambda)') +
  theme(legend.position = "none")

## Warning: Removed 5 rows containing missing values (geom_errorbar).

## Warning: Removed 4 rows containing missing values (geom_path).

(lasso_lowest_rmse <- lasso_grid %>%
  select_best("rmse", maximize = FALSE))

## # A tibble: 1 x 1
##   penalty
##     <dbl>
## 1 0.00359

Modelo Lasso final: O pacote tidymodels permite finalizar o modelo com o melhor parâmetro encontrado. Após isso, o modelo é aplicado nas amostras criadas com cross-validatione também ajustada na base de treino completa.

lasso_final <- lasso_spec %>% 
  finalize_model(parameters = lasso_lowest_rmse)

lasso_res <- fit_resamples(Sale_Price ~ .,
               lasso_final,
               cv_splits,
               control = control_resamples(save_pred = TRUE))
    
(lasso_res %>% 
  collect_metrics())

## # A tibble: 2 x 5
##   .metric .estimator  mean     n std_err
##   <chr>   <chr>      <dbl> <int>   <dbl>
## 1 rmse    standard   0.340     5  0.0289
## 2 rsq     standard   0.883     5  0.0171

lasso_fit <- lasso_final %>% 
      fit(Sale_Price ~.,
      data = train_baked)

Infelizmente essa forma de ajustar os modelos não permite a visualização das variáveis pela biblioteca plotmo, pois a classe criada pelo tidymodels não é reconhecida. É possível, no entanto, montar o gráfico, mas nesse caso específico, devido à quantidade de variáveis, a visualização fica prejudicada.

#library(plotmo)
#plot(lasso_fit, xvar = 'lambda', label = FALSE)

#?plot_glmnet

### Fazer fit na base inteira
#plot_glmnet(lasso_fit)

lasso_fit %>% 
  tidy %>% 
  mutate(log.lambda = log(lambda)) %>% 
  #filter(step <=100) %>% 
  #filter(estimate < 0.2) %>% 
  ggplot(aes(log.lambda, estimate, color = term)) +
  geom_line(show.legend = FALSE)

Será utilizado, portanto, o pacote vip, que possui funcionalidades semelhantes e permite identificar as variáveis mais importantes no modelo.

#plot(lasso.coef[lasso.coef != 0])
lasso_var <- lasso_fit %>%
  vi(lambda = lasso_lowest_rmse$penalty) %>% 
  mutate(Importance_pct = abs(Importance)/max(abs(Importance))) %>% 
  mutate(Variable = fct_reorder(Variable, Importance_pct))

#Verificaçãoda seleção das variáveis
lasso_var %>% 
  count(Importance_pct != 0)

## # A tibble: 2 x 2
##   `Importance_pct != 0`     n
##   <lgl>                 <int>
## 1 FALSE                     4
## 2 TRUE                    185

Aparentemente o modelo removeu apenas 4 variáveis, o que não ajuda muito nesse caso.

#variáveis que mais impactam no preço  
lasso_var %>% 
  filter(Importance_pct > 0.05) %>% 
  ggplot(aes(Variable,Importance_pct, fill = Sign)) +
  geom_col()+
  scale_y_continuous(labels = scales::percent_format())+
  coord_flip()

#variáveis que impactam negativamente no preço
#lasso_var %>% 
#  filter(Importance < 0) %>% 
#  ggplot(aes(Variable,Importance_pct, fill = Sign)) +
#  geom_col()+
#  scale_y_continuous(labels = scales::percent_format())+
#  coord_flip()

Observa-se que as principais variáveis são da mesma categoria e não fazem muito sentido se compararmos com o que foi visto na análise exploratória. Talvez rodar o modelo novamente com outros parâmetros e pré-processamento deva ser considerado. Foi filtrado o percntual acima de 5% para o gráfico não ficar muito poluído.

4.3. Ridge Regression

Assim, como no Lasso, o lambda será otimizado.

ridge_spec <- linear_reg(penalty = tune(), mixture = 0) %>% #mixture = 0 para ridge regression
  set_engine("glmnet", standardize = FALSE)

ridge_lambda_grid <- grid_regular(penalty(), levels = 1000)

doParallel::registerDoParallel() #processamento em paralelo para otimizar

set.seed(2020)
ridge_grid <- tune_grid(Sale_Price ~.,
  model = ridge_spec,
  resamples = ames_boot,
  grid = ridge_lambda_grid
)

ridge_grid %>%
  collect_metrics()

## # A tibble: 2,000 x 6
##     penalty .metric .estimator  mean     n std_err
##       <dbl> <chr>   <chr>      <dbl> <int>   <dbl>
##  1 1.00e-10 rmse    standard   0.355    10 0.0120 
##  2 1.00e-10 rsq     standard   0.879    10 0.00728
##  3 1.02e-10 rmse    standard   0.355    10 0.0120 
##  4 1.02e-10 rsq     standard   0.879    10 0.00728
##  5 1.05e-10 rmse    standard   0.355    10 0.0120 
##  6 1.05e-10 rsq     standard   0.879    10 0.00728
##  7 1.07e-10 rmse    standard   0.355    10 0.0120 
##  8 1.07e-10 rsq     standard   0.879    10 0.00728
##  9 1.10e-10 rmse    standard   0.355    10 0.0120 
## 10 1.10e-10 rsq     standard   0.879    10 0.00728
## # … with 1,990 more rows

ridge_grid %>%
  collect_metrics() %>%
  ggplot(aes(penalty, mean, color = .metric)) +
  geom_errorbar(aes(
    ymin = mean - std_err,
    ymax = mean + std_err
  ),
  alpha = 0.5
  ) +
  geom_line(size = 1.5) +
  facet_wrap(~.metric, scales = "free", nrow = 2) +
  scale_x_log10() +
  theme(legend.position = "none")

(ridge_lowest_rmse <- ridge_grid %>%
  select_best("rmse", maximize = FALSE))

## # A tibble: 1 x 1
##   penalty
##     <dbl>
## 1   0.117

Modelo Final:

ridge_final <- ridge_spec %>% 
  finalize_model(parameters = ridge_lowest_rmse)

ridge_res <- fit_resamples(Sale_Price ~ .,
               ridge_final,
               cv_splits,
               control = control_resamples(save_pred = TRUE))

(ridge_res %>% 
  collect_metrics())

## # A tibble: 2 x 5
##   .metric .estimator  mean     n std_err
##   <chr>   <chr>      <dbl> <int>   <dbl>
## 1 rmse    standard   0.340     5  0.0273
## 2 rsq     standard   0.884     5  0.0159

ridge_fit <- ridge_final %>% 
      fit(Sale_Price ~.,
      data = train_baked)

ridge_fit %>% 
  summary

##         Length Class      Mode   
## lvl      0     -none-     NULL   
## spec     6     linear_reg list   
## fit     12     elnet      list   
## preproc  5     -none-     list   
## elapsed  5     proc_time  numeric

ridge_fit %>% 
  tidy %>% 
  mutate(log.lambda = log(lambda)) %>% 
  #filter(step <=100) %>% 
  #filter(estimate < 0.2) %>% 
  ggplot(aes(log.lambda, estimate, color = term)) +
  geom_line(show.legend = FALSE)

#plot(lasso.coef[lasso.coef != 0])
var_ridge <- ridge_fit %>%
  vi(lambda = ridge_lowest_rmse$penalty) %>% 
  mutate(Importance_pct = abs(Importance)/max(abs(Importance))) %>% 
  mutate(Variable = fct_reorder(Variable, Importance_pct))

#variáveis que mais impactam no preço  
var_ridge %>% 
  filter(Importance_pct > 0.20) %>% 
  ggplot(aes(Variable,Importance_pct, fill = Sign)) +
  geom_col()+
  scale_y_continuous(labels = scales::percent_format())+
  coord_flip()

#variáveis que impactam negativamente no preço
#var_ridge %>% 
#  filter(Importance < 0) %>% 
#  ggplot(aes(Variable,Importance_pct, fill = Sign)) +
#  geom_col()+
#  scale_y_continuous(labels = scales::percent_format())+
#  coord_flip()

Apesar de não zerar o coeficiente de nenhuma variável, o modelo ridge selecionou varíveis que estão mais correlacionadas com a variável resposta. Além disso, observa-se o impacto negativo de algumas das principais selecionadas (mais de 50%).

4.4. Bagging

Assim como a seleção stepwise, também não há uma função específica no tidymodels para bagging. Seria possível através do crescimento de várias árvores, mas para simplificar, será utilizado o pacote ipred que realiza esse procedimento muito mais rápido.

library(ipred)

n_trees = seq(50,2200,500)

#função para calcular o eqm
f_eqm <- function(model, y = train_baked$Sale_Price){
          sum((y - round(model$y))^2) / length(model$y)
        }

bag_res <- tibble(trees = n_trees,
                  oob = NA,
                  rmse = NA)

for (i in 1:nrow(bag_res)){
  
  #print(i)
  
  set.seed(123)
  bag_fit <- ipred::bagging(Sale_Price ~ ., data = train_baked, coob = TRUE,  nbagg = n_trees[i])
  bag_res$oob[i] <- bag_fit$err
  bag_res$rmse[i] <- f_eqm(bag_fit)
}

A partir dos resultados computados, avalia-se o erro out-of-bag para diferentes valores de árvores:

bag_res %>% 
  ggplot(aes(trees,oob))+
  geom_line()

Como podemos observar, mesmo sendo uma diferença pequena, o erro estabiliza a partir de 1000 “bagged” árvores. Com isso, temos o modelo final.

bag_fit <- ipred::bagging(Sale_Price ~ ., data = train_baked, coob = TRUE, nbagg = 500)

4.5. Random Forest

Aqui será feita a otimização do parâmetro mtry e o número de árvores. Para otimizar o processamento, iniciaremos com um grid apenas para o mtry, para depois testar o número de árvores.

p <- ncol(train_baked) - 1 #total de variáveis preditoras

rf_tune_spec <- rand_forest(mode = "regression",
                            mtry = tune(), #p/3 = 189
                            trees = 500) %>% 
                set_engine("ranger", importance = "permutation")


(rf_grid <- grid_regular(
  mtry(range = c(p/5, 100)),
  #trees(range = c(500, 2000)),
  levels = 5))

## # A tibble: 5 x 1
##    mtry
##   <int>
## 1    37
## 2    53
## 3    68
## 4    84
## 5   100

doParallel::registerDoParallel() #processamento em paralelo para otimizar

set.seed(2020)
rf_tune <- tune_grid(Sale_Price ~.,
                      model = rf_tune_spec,
                      resamples = ames_boot,
                      grid = rf_grid)

Avaliando os resultados:

rf_tune %>% 
  collect_metrics()

## # A tibble: 10 x 6
##     mtry .metric .estimator  mean     n std_err
##    <int> <chr>   <chr>      <dbl> <int>   <dbl>
##  1    37 rmse    standard   0.326    10 0.0117 
##  2    37 rsq     standard   0.904    10 0.00630
##  3    53 rmse    standard   0.325    10 0.0114 
##  4    53 rsq     standard   0.903    10 0.00619
##  5    68 rmse    standard   0.325    10 0.0115 
##  6    68 rsq     standard   0.903    10 0.00632
##  7    84 rmse    standard   0.326    10 0.0115 
##  8    84 rsq     standard   0.901    10 0.00642
##  9   100 rmse    standard   0.327    10 0.0117 
## 10   100 rsq     standard   0.900    10 0.00662

rf_tune%>%
  collect_metrics() %>%
  select(mean, mtry, .metric) %>%
  filter(.metric == 'rmse') %>% 
  #pivot_longer(min_n:mtry,
  #  values_to = "value",
  #  names_to = "parameter"
  #) %>%
  ggplot(aes(mtry, mean, color = .metric)) +
  geom_point(show.legend = TRUE) #+

  #facet_wrap(~parameter, scales = "free_x") +
  #labs(x = NULL, y = "Value")

rf_lowest_rmse <- rf_tune %>%
  select_best("rmse", maximize = FALSE)

best_mtry = rf_lowest_rmse$mtry

Observa-se que mtry ótimo está perto de 50. Acima disso, o modelo provavelmente está overfitting. Com o melhor mtry, vamos avaliar o erro considerando o número de árvores e finalizar o modelo com os melhores parâmetros.

Melhor número de árvores

rf_spec <- rand_forest(mode = "regression",
                            mtry = best_mtry,
                            trees = tune()) %>% 
                set_engine("ranger", importance = "permutation")


rf_grid <- grid_regular(
  #mtry(range = c(p/4,100)),
  trees(range = c(500, 2000)),
  levels = 4)

doParallel::registerDoParallel()

set.seed(2020)
rf_tune <- tune_grid(Sale_Price ~.,
                      model = rf_spec,
                      resamples = ames_boot,
                      grid = rf_grid)

rf_tune %>% 
  collect_metrics()

## # A tibble: 8 x 6
##   trees .metric .estimator  mean     n std_err
##   <int> <chr>   <chr>      <dbl> <int>   <dbl>
## 1   500 rmse    standard   0.325    10 0.0115 
## 2   500 rsq     standard   0.903    10 0.00629
## 3  1000 rmse    standard   0.324    10 0.0114 
## 4  1000 rsq     standard   0.904    10 0.00621
## 5  1500 rmse    standard   0.324    10 0.0113 
## 6  1500 rsq     standard   0.904    10 0.00624
## 7  2000 rmse    standard   0.324    10 0.0114 
## 8  2000 rsq     standard   0.904    10 0.00625

rf_tune %>%
  collect_metrics() %>%
  select(mean, trees, .metric) %>%
  filter(.metric == 'rmse') %>% 
  #pivot_longer(min_n:mtry,
  #  values_to = "value",
  #  names_to = "parameter") %>%
  ggplot(aes(trees, mean, color = .metric)) +
  geom_point(show.legend = TRUE)

  rf_lowest_rmse <- rf_tune %>%
    select_best("rmse", maximize = FALSE)

Modelo Final:

rf_final <-  rf_spec %>%
    finalize_model(parameters = rf_lowest_rmse)


rf_res <- fit_resamples(Sale_Price ~ .,
               rf_final,
               cv_splits,
               control = control_resamples(save_pred = TRUE))


rf_fit <- rf_final %>% 
      fit(Sale_Price ~.,
      data = train_baked)

rf_fit %>% 
  summary

##         Length Class       Mode   
## lvl      0     -none-      NULL   
## spec     6     rand_forest list   
## fit     15     ranger      list   
## preproc  1     -none-      list   
## elapsed  5     proc_time   numeric

Importância das variáveis

vi(rf_fit) %>% 
  mutate(Importance_pct = abs(Importance)/max(abs(Importance))) %>% 
  mutate(Variable = fct_reorder(Variable, Importance_pct)) %>% 
  filter(Importance_pct > 0.05) %>% 
  ggplot(aes(Variable, Importance_pct)) +
  geom_point()+
  scale_y_continuous(labels = scales::percent_format())+
  coord_flip()

5. Comparação dos Modelos

No gráfico abaixo, é possível observar o erro quadrado médio (rmse) e o R-quadrado (rsq) para cada um dos modelos, exceto para bagging, que foi ajustado diretamente na base completa.

lm_res %>% 
  select(id, .metrics) %>% 
  unnest(.metrics) %>% 
  mutate(model = "linear regression") %>% 
  bind_rows(lasso_res %>% 
            select(id, .metrics) %>% 
            unnest(.metrics) %>% 
            mutate(model = "lasso")) %>% 
   bind_rows(ridge_res %>% 
            select(id, .metrics) %>% 
            unnest(.metrics) %>% 
            mutate(model = "ridge")) %>% 
  #bind_rows(bag_res %>% 
   #         select(id, .metrics) %>% 
    #        unnest(.metrics) %>% 
     #       mutate(model = "bagging")) %>% 
  bind_rows(rf_res %>% 
            select(id, .metrics) %>% 
            unnest(.metrics) %>% 
            mutate(model = "random forest")) %>% 
  ggplot(aes(id, .estimate, group = model, color = model)) + 
    geom_point(size = 1.5) + 
    facet_wrap(~.metric) + 
    coord_flip()

Para ter uma ideia melhor, segue a comparação das métricas nas bases de treino e teste completas.

results_train <- lm_fit %>% 
  predict(new_data = train_baked) %>% 
  mutate(truth = train_baked$Sale_Price,
         model = 'lm') %>% 
  bind_rows(lasso_fit %>% 
  predict(new_data = train_baked) %>% 
  mutate(truth = train_baked$Sale_Price,
         model = 'lasso')) %>%
  bind_rows(ridge_fit %>% 
  predict(new_data = train_baked) %>% 
  mutate(truth = train_baked$Sale_Price,
         model = 'ridge')) %>%
  bind_rows(tibble(.pred = predict(bag_fit, newdata = train_baked),
                    truth = train_baked$Sale_Price,
                    model = 'bagging')) %>% 
  bind_rows(rf_fit %>% 
  predict(new_data = train_baked) %>% 
  mutate(truth = train_baked$Sale_Price,
         model = 'random forest'))

results_train %>% 
  group_by(model) %>% 
  rmse(truth = truth, estimate = .pred)

## # A tibble: 5 x 4
##   model         .metric .estimator .estimate
##   <chr>         <chr>   <chr>          <dbl>
## 1 bagging       rmse    standard       0.393
## 2 lasso         rmse    standard       0.290
## 3 lm            rmse    standard       0.270
## 4 random forest rmse    standard       0.128
## 5 ridge         rmse    standard       0.295

Pela base de treino, observa-se que o melhor modelo foi random forest, seguido pela regressão linear. Para confirmar, vamos observar o comportamento na base teste:

results_test <- lm_fit %>% 
  predict(new_data = test_baked) %>% 
  mutate(truth = test_baked$Sale_Price,
         model = 'lm') %>% 
  bind_rows(lasso_fit %>% 
  predict(new_data = test_baked) %>% 
  mutate(truth = test_baked$Sale_Price,
         model = 'lasso')) %>%
  bind_rows(ridge_fit %>% 
  predict(new_data = test_baked) %>% 
  mutate(truth = test_baked$Sale_Price,
         model = 'ridge')) %>%
  bind_rows(tibble(.pred = predict(bag_fit, newdata = test_baked),
                    truth = test_baked$Sale_Price,
                    model = 'bagging')) %>% 
  bind_rows(rf_fit %>% 
  predict(new_data = test_baked) %>% 
  mutate(truth = test_baked$Sale_Price,
         model = 'random forest'))

## Warning in predict.lm(object = object$fit, newdata = new_data, type =
## "response"): prediction from a rank-deficient fit may be misleading

results_test %>% 
  group_by(model) %>% 
  rmse(truth = truth, estimate = .pred)

## # A tibble: 5 x 4
##   model         .metric .estimator .estimate
##   <chr>         <chr>   <chr>          <dbl>
## 1 bagging       rmse    standard       0.481
## 2 lasso         rmse    standard       0.503
## 3 lm            rmse    standard       0.557
## 4 random forest rmse    standard       0.375
## 5 ridge         rmse    standard       0.479

Na base de teste, tivemos o mesmo resultado: o modelo que melhor performou foi random forest. No entanto, o segundo melhor modelo foi bagging, que teve a pior performance na base de treino. Ridge e Lasso também tiveram a performance invertida. A regressão linear, no entanto, foi o modelo que performou pior nesse caso.

6. Conclusão

A floresta aleatória foi o modelo que melhor performou, tanto na base treino, como na base teste e, portanto, o modelo selecionado como melhor. Conforme avaliado na análise exploratória, as principais variáveis preditoras estavam altammente correlacionadas, mas não necessariamente no topo da lista. As 3 principais foram: Overall_Quality, Gr_Liv_Area e Year_Built.

7. Referências

Como referências para o desenvolvimento do relatório foram utilizadas as seguintes fontes:

An Introduction do Statistical Modelling - James, G. et. al.
Material de aula do curso Modelos Preditivos - Insper - 2020
Introdução a Tidymodels - Mendonça, T.
Julia Silge - Blog e YouTube

Bootstrap resampling with #TidyTuesday beer production data

Preprocessing and resampling using #TidyTuesday college data

LASSO regression using tidymodels and #TidyTuesday data for The Office

Tuning random forest hyperparameters with #TidyTuesday trees data

#TidyTuesday hotel bookings and recipes

#TidyTuesday and tidymodels
AmesHousing
Tidymodels documentation
Kaggle: House Prices: Glmnet, XGBoost, and SVM Using tidymodels House prices: Lasso, XGBoost, and a detailed EDA
Hands on Machine Learning with R

AmesHousing com tidymodels

Viviane Sanchez

4/9/2020

1. Introdução

Objetivo:

Bibliotecas

Dados

2. Análise Exploratória

2.1. Distribuição dos preços

2.2. Correlação entre variáveis numéricas

2.3. Ano de construção, Área total e Material de construção

2.4. Garagem

2.5. Área total, Lareiras e Condição de Venda

2.5 Qualidade de venda

3. Pré processamento

3.1 Separação da base em treino e teste (treino/teste/total):

3.2. Receita

4. Modelos

4.1. Regressão Linear

4.1.1. Sem seleção stepwise

4.1.2 Com seleção stepwise backward

4.2. Lasso

4.3. Ridge Regression

4.4. Bagging

4.5. Random Forest

5. Comparação dos Modelos

6. Conclusão

7. Referências